{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bayesian Data Analysis, 3rd ed\n",
    "## Chapter 2, demo 3\n",
    "\n",
    "Authors:\n",
    "- Aki Vehtari <aki.vehtari@aalto.fi>\n",
    "- Tuomas Sivula <tuomas.sivula@aalto.fi>\n",
    "\n",
    "Probability of a girl birth given placenta previa (BDA3 p. 37).\n",
    "Simulate samples from Beta(438,544), draw a histogram with quantiles, and do the same for a transformed variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# import necessary packages\n",
    "\n",
    "import numpy as np\n",
    "from scipy.stats import beta\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# add utilities directory to path\n",
    "import os, sys\n",
    "util_path = os.path.abspath(os.path.join(os.path.pardir, 'utilities_and_data'))\n",
    "if util_path not in sys.path and os.path.exists(util_path):\n",
    "    sys.path.insert(0, util_path)\n",
    "\n",
    "# import from utilities\n",
    "import plot_tools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# edit default plot settings\n",
    "plt.rc('font', size=12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAH1CAYAAACtCxSGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHWdJREFUeJzt3X20bWVdL/DvDw4vgwxILMwQQTuIL4AOvd0KEwVCyGvm\nS6UxFM2ujbymRGC3fIlK8ZaVpqZdGeZLmuQwDYOkMiSVvAloceEmYICBhEoWAr6k8Nw/5jqDxT77\nvOz9rLXn2vt8PmOscfaaa665f/M5e67vfOZ81pzVWgsAsHq7jV0AAKx3whQAOglTAOgkTAGgkzAF\ngE7CFAA6CVMA6CRMN4Cq2quq3lpVn6uq26rqH6rqpG3M+5yqurOqbp96PG7y2qaqOqeq/qOqLqiq\nfafe9ytVddoarRIshKp6SFVdWFW3VtVnq+opk+knL9mGvlpVraoetY3lXFRVX5+a/6qp146qqiur\n6pbpbayq9qiqv6+q+89/TeklTDeGTUluSHJMkv2SvCzJe6vqkG3M/4nW2r2mHhdNpj81SUtynyS3\nJnl+klTVoUl+NMnr57UCsGiqalOSc5Ocl+TeGbaHd1XVYa21d09vQ0lekOTaJJ/aziJfOPWeB09N\nf3WS05McleSlVXXfyfTTkvxpa+2GGa8acyBMN4DW2h2ttTNba9e31u5qrZ2X5Loky+4lb8ehSS5q\nrX0ryUeSPHAy/fVJfnEyHXYVhye5X5LXttbubK1dmOTiJM9aZt5Tkryzre6ScocmubC19vkk1yQ5\nuKoekORpSV67utJZa8J0A6qqA5McluTKbczyyMkhpaur6uWTPfAkuSLJsVW1V5LHJ7lycljrltba\nxfOvHBZeJXn4PSYMwffYJO/cwXtfPdnuLt5yamXiiiQnVNVBSQ5J8s9Jfi/JGa21b86qcOZLmG4w\nVbVHkncneUdr7TPLzPLRDB8G35Vhz/eZSc6YvPYXGXq0l2Q4zHtOkl9N8pKqelVVfbSq3lRVe855\nNWARXJXki0nOmJy/PCHDqZR9lsz37CQfa61dt51l/VKGIz3fk+QtSf68qh40ee30JD+X5INJfiHJ\n0UluS3JdVZ1bVX9bVT8+q5ViPsqF7jeOqtotyR8n2TfJk3dmr7aqnpFhD3irQ8JV9ZoMe8n/kuTF\nSU5McnaSS1trfzDL2mERVdWRSd6QYQf00iRfSvKN1trzpua5JslZrbW3rWC5FyQ5v7X2hiXT90ny\niSQnTH7vnyU5P0Pv9ajW2pf71oh50TPdIKqqkrw1yYFJnraCw0Mtw6Grpcs7IskPZtiLPiLJZZPz\nQZckOXImRcOCa61d3lo7prV2QGvtCRl6l5/c8npVHZ3hvOr7VrroLLPdJXlFkrNba1/IsN1d2lq7\nNcmNSb53NevA2hCmG8ebkzwkyZNaa1/b1kxVddLknGqq6vAkL88wYnF6nkryxiQvaq3dleHQ72Mm\nh3ePyTBqETa8qjqyqvauqn2q6vQk353k7VOznJJhxO1t21nG/lX1hMlyNlXVyRnOsV6wZL6HJnlc\nhm05Gba7Yyfb6+YMR4hYUMJ0A5gMgPjZJI9IcvPUd9lOrqqDJz8fPJn9uCSXV9UdGc6Rvj/JWUsW\n+dwkV7TWLps8f3+SmzIc4jogQ28VdgXPSvKvGc6dHpfkh1tr30iSqto7yU8kecfSN02+l/2hydM9\nkrwyw/ZzS5KfT/JjrbWrl7zt95O8uLV25+T5Lyd5UYaBhGe11m6e5YoxW86ZAkAnPVMA6LRpx7Pc\ng27sRvG2Jw7/Pvf8ceuAjcx2thEsN1BsK3qmANBJmAJAJ2EKAJ2EKQB0EqYA0EmYAkAnYQoAnVb6\nPVNgBK/78DUrfs+px2+eQyXAcvRMAaCTMAWATsIUADoJUwDoJEwBoJMwBYBOwhQAOglTAOgkTAGg\nkzAFgE4uJwhraDWXBQQWn54pAHQSpgDQSZgCQCdhCgCdhCkAdDKaF1bJyFxgCz1TAOgkTAGgkzAF\ngE7CFAA6CVMA6GQ0L2xQqx1tfOrxm2dcCWx8eqYA0EmYAkAnYQoAnYQpAHQyAAm4h9UMXDJoiV2d\nnikAdBKmANBJmAJAJ2EKAJ2EKQB0EqYA0EmYAkAnYQoAnYQpAHRyBSR2eau9VRnAFnqmANBJmAJA\nJ2EKAJ2EKQB0EqYA0EmYAkAnYQoAnYQpAHQSpgDQSZgCQCdhCgCdhCkAdHKhe6Dbam8WcOrxm2dc\nCYxDzxQAOglTAOjkMC8binuTAmPQMwWATsIUADoJUwDoJEwBoJMwBYBOwhQAOglTAOgkTAGgkzAF\ngE7CFAA6CVMA6CRMAaCTMAWATsIUADoJUwDoJEwBoJMwBYBOwhQAOglTAOgkTAGgkzAFgE7CFAA6\nCVMA6CRMAaDTprELAHZdr/vwNat636nHb55xJdBHzxQAOglTAOgkTAGgkzAFgE7CFAA6CVMA6CRM\nAaCTMAWATsIUADq5AhILabVXxgEYg54pAHQSpgDQSZgCQCdhCgCdhCkAdBKmANBJmAJAJ2EKAJ2E\nKQB0cgUkYN1ZzRWyTj1+8xwqgYGeKQB0EqYA0EmYAkAnYQoAnYQpAHQSpgDQSZgCQCdhCgCdhCkA\ndBKmANDJ5QSZu9Vc+g1gPdEzBYBOwhQAOglTAOgkTAGgkzAFgE7CFAA6CVMA6OR7psAuYbXfdz71\n+M0zroSNSM8UADoJUwDoJEwBoJMwBYBOwhQAOhnNy05z9xeA5emZAkAnYQoAnYQpAHQSpgDQSZgC\nQCdhCgCdhCkAdPI9013Ujf/+1STJ+3x3FKCbnikAdBKmANBJmAJAJ2EKAJ2EKQB0MpoXYDtWc7ek\nU4/fPIdKWGR6pgDQSZgCQCdhCgCdhCkAdBKmANDJaF6AGdsyAvjpK7wGtlHA65eeKQB0EqYA0EmY\nAkAn50zXudVcnSVJnj7jOgB2ZXqmANBJmAJAJ2EKAJ2EKQB0MgAJYEG43dv6pWcKAJ2EKQB0EqYA\n0Mk50wWy2gswADAuPVMA6KRnOgd6mAC7FmEKsI6tdufdV2pmy2FeAOgkTAGgk8O8O+D8J7ARudrS\nbOmZAkCnXaZnqocJ0Gc9fI6O1Xuu1trOz1x1QZL7zK+cdec+SW4Zu4gFoj22pk3uSXtsTZtsbZHa\n5JbW2ok7mmlFYco9VdWlrbVHj13HotAeW9Mm96Q9tqZNtrYe28Q5UwDoJEwBoJMw7fOWsQtYMNpj\na9rknrTH1rTJ1tZdmzhnCgCd9EwBoJMwBYBOwhQAOgnTKVV176r6QFXdUVWfq6qf2sH8e1bVP1XV\njUumP6KqLquqr07+fcR8K5+fWbRJVR1WVedW1Zeq6stV9ZdV9eD5Vz97s/obmXr92VXVqupn5lPx\n/M1wu9m9ql5ZVTdV1W1V9emq2n++1c/HDNvk2Kr6VFV9paqurarnz7fy+djZ9qiqM6vqm1V1+9Tj\ngVOvL+xnqzC9p99P8p9JDkxycpI3V9XDtjP/GUm+ND2hqvZMcm6SdyX5jiTvSHLuZPp61N0mSfZP\n8sEkD54s55MZ2mg9mkV7JEmq6juS/EqSK2dd5BqbVZv8WpIfTPIDSfZN8qwkX59tqWtmFp8leyT5\nQJL/nWS/JD+Z5Her6qi5VDxfK2mPP2mt3WvqcW2yDj5bW2sew4jmb8vwn33Y1LQ/SvK/tjH/oUn+\nKclJSW6cmn5Cks9nMlJ6Mu1fkpw49jqO1SbLzHfvJC3JAWOv45jtkeQPkrwgyUVJfmbs9RuzTTJ8\nON6e5EFjr9MCtcmBk+1kn6lplyR55tjrOK/2SHJmkndtYzkL/dmqZ3q3w5J8q7V29dS0f0yyrb2n\nN2ToVXxtyfSHJbm8Tf6nJy7fznIW2azaZKnHJrm5tfZv/SWuqZm1R1V9X5JHZwjU9WxWbXJEkm8l\neXpV3VxVV1fV/5h5tWtjJm3SWvtCkvckee7kEPgPJHlAko/PvuS5Wml7PGlyOujKqvq5qekL/dkq\nTO92ryRfWTLt1iTfvnTGqnpKkt1bax/YxnJu3ZnlrAOzapPp+Q7KcMjntFkVuYZm0h5VtXuSNyV5\nYWvtrnkUuoZm9TdyUIZDmYdl6Kk9PcmZVfXDsy13Tcxyu3lPklck+UaSjyV5aWvthhnWuhZ2uj2S\nvDfJQ5J8Z5L/nuQVVfXMqeUs7GfrLnMLtp1we4bzNNP2TXLb9ISq+rYkv5XkR3qWs07Mqk22zPed\nSf4qyZtaa++ZYZ1rZVbt8YIMe9j/Z+YVrr1ZtcmWXtmvt9a+luTyqjpnMv9fz67cNTGTNqmqw5Oc\nk+SpGdpgc5Lzquqm1tr5sy56jnb6M7G19v+mnv5dVf1ehh2r96xkOWPQM73b1Uk2VdX0zfCOytaD\nQzYnOSTJx6rq5iTvT/Ldk0NTh0zmP7Kqauo9Ry6znPVgVm2yZbDNXyX5YGvtVXOue15m1R7HJXnK\n5PnNGQbd/E5VvXHO9c/DrNrk8sl804fw1uvl2WbVJg9PcnVr7S9ba3e11q5Kcn6Gc6vryc62x3Ja\nki2fpYv92Tr2SdtFemTYC3xPhhPmR2c4hPCwJfNsSnLfqcdTk9w0+Xn3JHsm+VySFyfZK8kLJ8/3\nHHv9RmyTfTOM4H3j2OuzIO2x/5LX/y7DYe/9xl6/sdpkMs9HM4xc3SvDob4vJjlu7PUb8e/kQRl6\nY8dmCJQHJflskuePvX7zaI/JfE/OMBitknxfhgFHp0xeW+jP1tELWKRHhlGmf5bkjgyjxH5qMv2H\nkty+jfc8LktGaiZ5ZJLLMhy6+lSSR469bmO2SZJTMuxh3jH5cNjyOHjs9Rvrb2TJ6xdlnY7mnWWb\nJPmeJBdM/jauTfKzY6/bArTJTyS5IsOhzBuT/GaS3cZev3m1xyRw/23yN/CZJC9aspyF/Wx1oXsA\n6OScKQB0EqYA0EmYAkAnYQoAnYQpAHQSpgDQSZgCQCdhCgusqh5ZVRdPbob8yao6eOyagK0JU1hQ\nkzvs/EWGq94ckOGqQC8btShgWcIUFtfvJDm7tfbBNtxJ5Zwk/2XkmoBluAUbLKCq2jfDRb8Pm5q8\nW5Kvj1MRsD3CFBbTcUn2yHBfzy3T9kpy7mgVAdvkMC8spkMy3Pt1/y2PJB/JcFcVYMEIU1hMeyX5\n6pYnVXVokkcn+eBoFQHbJExhMV2S5Jiqul9V3T/JHyd5aWvtyyPXBSzDOVNYTBcmOS/J1Rlulvyb\nrbWzxy0J2BY3BweATg7zAkAnYQoAnYQpAHQSpgDQSZgCQCdhCgCdhCkAdBKmANBJmAJAJ2EKAJ2E\nKQB0EqYA0EmYAkAnYQoAnYQpAHQSpgDQSZgCQCdhCgCdhCkAdBKmANBJmAJAJ2EKAJ2EKQB0EqYA\n0EmYAkAnYQoAnYQpAHQSpgDQSZgCQCdhCgCdhCkAdBKmANBJmAJAJ2EKAJ2EKQB0EqYA0EmYAkAn\nYQoAnYQpAHQSpgDQSZiuoaraq6reWlWfq6rbquofquqkbcz7nKq6s6pun3o8bvLapqo6p6r+o6ou\nqKp9p973K1V12hqtEoyiqh5SVRdW1a1V9dmqespk+slLtpmvVlWrqkdtYzkXVdXXp+a/auq1o6rq\nyqq6ZXqbqqo9qurvq+r+819T1gthurY2JbkhyTFJ9kvysiTvrapDtjH/J1pr95p6XDSZ/tQkLcl9\nktya5PlJUlWHJvnRJK+f1wrA2KpqU5Jzk5yX5N4Z/v7fVVWHtdbePb3NJHlBkmuTfGo7i3zh1Hse\nPDX91UlOT3JUkpdW1X0n009L8qettRtmvGqsY8J0DbXW7mitndlau761dldr7bwk1yVZdq95Ow5N\nclFr7VtJPpLkgZPpr0/yi5PpsFEdnuR+SV7bWruztXZhkouTPGuZeU9J8s7WWlvF7zk0yYWttc8n\nuSbJwVX1gCRPS/La1ZXORiVMR1RVByY5LMmV25jlkZNDTFdX1csne+RJckWSY6tqrySPT3Ll5DDX\nLa21i+dfOSycSvLwe0wYgu+xSd65g/e+erKdXbzlVMrEFUlOqKqDkhyS5J+T/F6SM1pr35xV4WwM\nwnQkVbVHkncneUdr7TPLzPLRDB8O35VhT/iZSc6YvPYXGXq0l2Q4zHtOkl9N8pKqelVVfbSq3lRV\ne855NWAMVyX5YpIzJucvT8hw6mSfJfM9O8nHWmvXbWdZv5ThyM73JHlLkj+vqgdNXjs9yc8l+WCS\nX0hydJLbklxXVedW1d9W1Y/PaqVY32p1Rz/oUVW7JfnjJPsmefLO7OVW1TMy7BFvdUi4ql6TYa/5\nX5K8OMmJSc5Ocmlr7Q9mWTssgqo6MskbMuxwXprkS0m+0Vp73tQ81yQ5q7X2thUs94Ik57fW3rBk\n+j5JPpHkhMnv/bMk52fovR7VWvty3xqx3umZrrGqqiRvTXJgkqet4HBRy3Aoa+nyjkjygxn2qo9I\nctnk/NAlSY6cSdGwYFprl7fWjmmtHdBae0KG3uUnt7xeVUdnOK/6vpUuOstsZ0lekeTs1toXMmxn\nl7bWbk1yY5LvXc06sLEI07X35iQPSfKk1trXtjVTVZ00Oaeaqjo8ycszjGCcnqeSvDHJi1prd2U4\n9PuYyeHdYzKMYoQNp6qOrKq9q2qfqjo9yXcnefvULKdkGHF723aWsX9VPWGynE1VdXKGc6wXLJnv\noUkel2HbTYbt7NjJ9rk5wxEhdnHCdA1NBkT8bJJHJLl56rttJ1fVwZOfD57MflySy6vqjgznSN+f\n5Kwli3xukitaa5dNnr8/yU0ZDnkdkKG3ChvRs5L8a4Zzp8cl+eHW2jeSpKr2TvITSd6x9E2T72F/\naPJ0jySvzLC93JLk55P8WGvt6iVv+/0kL26t3Tl5/stJXpRh4OBZrbWbZ7lirE/OmQJAJz1TAOi0\nacez3INu7Fp72xOHf597/rh1wEZj22LnLDcgbSt6pgDQSZgCQCdhCgCdhCkAdBKmANBJmAJAJ2EK\nAJ1W+j1T2HBe9+FrVvW+U4/fPONKgPVKzxQAOglTAOgkTAGgkzAFgE7CFAA6CVMA6CRMAaCTMAWA\nTsIUADoJUwDoJEwBoJMwBYBOLnQPq7SaC+S7OD5sTHqmANBJmAJAJ2EKAJ2EKQB0EqYA0MloXjaU\n1YywBeilZwoAnfRMYQ2ttufs+6mw2PRMAaCTnikLyblPYD3RMwWATsIUADoJUwDoJEwBoJMwBYBO\nwhQAOglTAOgkTAGgkzAFgE7CFAA6CVMA6CRMAaCTMAWATsIUADoJUwDoJEwBoJMwBYBOwhQAOglT\nAOi0aewCgB173YevWfF7Tj1+8xwqAZajZwoAnYQpAHQSpgDQSZgCQCdhCgCdhCkAdBKmANDJ90yZ\nu9V8RxJgPdEzBYBOwhQAOglTAOgkTAGgkzAFgE7CFAA6CVMA6CRMAaCTizbABrXai2W4qTisnJ4p\nAHQSpgDQSZgCQCdhCgCdhCkAdBKmANBJmAJAJ2EKAJ2EKQB0EqYA0EmYAkAnYQoAnYQpAHQSpgDQ\nSZgCQCdhCgCdhCkAdBKmANBp09gFsH687sPXjF0CwELSMwWATsIUADoJUwDoJEwBoJMwBYBOwhQA\nOglTAOgkTAGgkzAFgE7CFAA6CVMA6CRMAaCTMAWATu4aA9zDau4OdOrxm+dQCawfeqYA0EmYAkAn\nYQoAnYQpAHQSpgDQSZgCQCdhCgCdfM90F7Wa7xICsDw9UwDoJEwBoJMwBYBOwhQAOglTAOgkTAGg\nkzAFgE7CFAA6CVMA6CRMAaCTMAWATq7NC3Rb7bWeTz1+84wrgXHomQJAJ2EKAJ2EKQB0EqYA0EmY\nAkAnYQoAnYQpAHQSpgDQSZgCQCdhCgCdhCkAdHJt3nVutddEBWB29EwBoJOeKTAad5tho9AzBYBO\nwhQAOglTAOgkTAGgkzAFgE7CFAA6+WoMsO6s5is1vk7DPOmZAkAnYQoAnYQpAHQSpgDQyQCkBbLc\noIqn//tXkyTvc3cYgIWlZwoAnYQpAHQSpgDQSZgCQCcDkIBdwtIBfjs7uM+Vk9gZeqYA0EmYAkAn\nYQoAnYQpAHQSpgDQyWjeOVjNvRYBWL/0TAGgkzAFgE7CFAA6OWcKsB2rGQPhqkm7Hj1TAOgkTAGg\nkzAFgE7CFAA6GYC0Ay7AAMCOCFOAGVvtTrhRwOuXMAVYEL6Gs345ZwoAnXaZnqlznwDMyy4TpgAb\nkfOzi2HUMNVbBBjHon/+rrewr9bazs9cdUGS+8yvnB26T5JbRvz9i0RbDLTD3bTF3bTF3bTFYLXt\ncEtr7cQdzbSiMB1bVV3aWnv02HUsAm0x0A530xZ30xZ30xaDebeD0bwA0EmYAkCn9Rambxm7gAWi\nLQba4W7a4m7a4m7aYjDXdlhX50wBYBGtt54pACwcYQoAnYQpAHRaqDCtqhdW1aVV9Y2qevsO5v2F\nqrq5qr5SVX9YVXutUZlzt7PtUFUPr6q/rKpbqmpDnvxeQVucUlWXTf4ebqyq36qqDXW5zBW0xTOq\n6qqqurWqvlhV76iqfdew1LlbyWfF1Hv+pqraLvx38ZyqurOqbp96PG7tKp2vFebHA6vqvKq6bfL5\n+Vu9v3+hwjTJTUlemeQPtzdTVT0hyf9MclySByR5YJJfm3t1a2en2iHJN5O8N8nz5l7ReHa2LfZJ\ncmqGq5z81wx/G6fPt7Q1t7NtcXGSo1tr+2XYNjZN3reR7GxbJEmq6uQke8y1ovGspC0+0Vq719Tj\novmWtqZ2Nj/2TPLXSS5Mct8kByV5V+8vX6g9tNba+5Okqh6dYQW35ZQkb22tXTmZ/zeSvDtDwK57\nO9sOrbWrklxVVd+7VrWttRW0xZunnn6+qt6d5PFzLm9NraAtblgy6c4kG+pvZAWfFamq/ZL8apJn\nJ/nE/KtbWytpi41sBe3wnCQ3tdZ+d2ra5b2/f9F6pjvrYUn+cer5PyY5sKoOGKkeFs9jk1w5dhFj\nqarHVNWtSW5L8rQkrxu5pDGdleTNSW4eu5AF8MjJYc2rq+rlG+2Q9076/iTXV9WHJm1xUVUd0bvQ\n9Rqm90py69TzLT9/+wi1sGCq6qeTPDrJb49dy1haax+fHOY9KMlrklw/bkXjmPRSjk7yhrFrWQAf\nTfLwJN+VYQfrmUnOGLWicRyU5BlJXp/kfknOT3Lu5PDvqq3XML09yfSAii0/3zZCLSyQqvqxJK9O\nclJrbZe/U0Zr7fNJLkhyzti1rLWq2i3Jm5K8uLX2rbHrGVtr7drW2nWttbtaa/83ya8nefrYdY3g\na0k+3lr7UGvtPzPsdB+Q5CE9C12vYXplkqOmnh+V5AuttX8bqR4WQFWdmOTsJE+afFgw2JTkQWMX\nMYJ9Mxyh+JOqujnJJZPpN1bVD41X1sJoSWrsIkZweYZ1n6mFCtOq2lRVeyfZPcnuVbX3No7pvzPJ\n86rqoVW1f5KXJXn7GpY6VzvbDjXYO8mek+d7b6SvCCUraotjMwxCe1pr7ZNrXedaWEFbnFxVB09+\nfkCSVyX5m7Wtdr52si1uzXAY7xGTx49Mpj8qyd+vWbFztoK/i5Oq6sDJz4cneXmSc9e22vlZQX68\nK8n3V9XxVbV7hm8B3JLkn7oKaK0tzCPJmRn2GKYfZyY5OMOh3YOn5j0tyReSfCXJ25LsNXb9a90O\nSQ5ZZr7rx65/pLb4SJJvTaZteXxo7PpHaotXJbkxyR2Tf9+S5ICx6x+jLZa8Z8v2smns+kf6u/jt\nyWfmHUmuzXCYd4+x6x/jbyLJU5N8dpIfFyV5WO/vd6F7AOi0UId5AWA9EqYA0EmYAkAnYQoAnYQp\nAHQSpgDQSZgCQCdhCgCdhCmsA1X1pqp6xdh1AMsTprA+PCbJx8cuAlieMIUFVVW7VdVLqur6JEck\nOaeqThu5LGAZwhQW1y8neWKSFyW5Osl/S/Kaqrr/qFUBWxGmsICq6tuTvDTJKUnun+TTbbi13A1J\nDh+zNmBrwhQW07FJrmqtXZ/kqCSfrqrdknxHki+OWRiwNWEKi+m+Sb48+fkRST6d5LEZbnh9+VhF\nActb7i7kwPg+k+RRVfXAJA/PEKxvT/KS5ibEsHDcHBwWUFVVkt9O8tNJ9ktyTZJXttb+aNTCgGUJ\nU1hgVfWTSZ7bWjtx7FqAbXPOFBbbg5NcNXYRwPYJU1hsD87wHVNggTnMCwCd9EwBoJMwBYBOwhQA\nOglTAOgkTAGgkzAFgE7CFAA6/X8HZJtJD4hPhwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa4e634bdd8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plotting grid\n",
    "x = np.linspace(0.36, 0.54, 150)\n",
    "\n",
    "# draw n random samples from Beta(438,544)\n",
    "n = 10000\n",
    "th = beta.rvs(438, 544, size=n)  # rvs comes from `random variates`\n",
    "\n",
    "# plot 2 subplots\n",
    "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(8, 8))\n",
    "# show only x-axis\n",
    "plot_tools.modify_axes.only_x(axes)\n",
    "# manually adjust spacing\n",
    "fig.subplots_adjust(hspace=0.5)\n",
    "\n",
    "# plot histogram\n",
    "axes[0].hist(th, bins=30, color=plot_tools.lighten('C0'))\n",
    "# compute 2.5% and 97.5% quantile approximation using samples\n",
    "th25, th975 = np.percentile(th, [2.5, 97.5])\n",
    "# draw lines for these\n",
    "axes[0].axvline(th25, color='C1')\n",
    "axes[0].axvline(th975, color='C1')\n",
    "axes[0].text(\n",
    "    th25,\n",
    "    axes[0].get_ylim()[1]+15,\n",
    "    '2.5%',\n",
    "    horizontalalignment='center'\n",
    ")\n",
    "axes[0].text(\n",
    "    th975,\n",
    "    axes[0].get_ylim()[1]+15,\n",
    "    '97.5%',\n",
    "    horizontalalignment='center'\n",
    ")\n",
    "axes[0].set_xlabel(r'$\\theta$')\n",
    "\n",
    "# plot histogram for the transformed variable\n",
    "phi = (1-th)/th\n",
    "axes[1].hist(phi, bins=30, color=plot_tools.lighten('C0'))\n",
    "# compute 2.5% and 97.5% quantile approximation using samples\n",
    "phi25, phi975 = np.percentile(phi, [2.5, 97.5])\n",
    "# draw lines for these\n",
    "axes[1].axvline(phi25, color='C1')\n",
    "axes[1].axvline(phi975, color='C1')\n",
    "axes[1].text(\n",
    "    phi25,\n",
    "    axes[1].get_ylim()[1]+15,\n",
    "    '2.5%',\n",
    "     horizontalalignment='center'\n",
    ")\n",
    "axes[1].text(\n",
    "    phi975,\n",
    "    axes[1].get_ylim()[1]+15,\n",
    "    '97.5%',\n",
    "    horizontalalignment='center'\n",
    ")\n",
    "axes[1].set_xlabel(r'$\\phi$');"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
